Ab-initio theory of metal-insulator interfaces in a finite electric field 
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We present a novel technique for calculating the dielectric response of metal/insulator heterostruc- 
tures. This scheme allows, for the first time, the fully first-principles calculation of the microscopic 
properties of thin-film capacitors at finite bias potential. The method can be readily applied to 
pure insulators, where it provides an interesting alternative to conventional finite-field techniques 
based on the Berry-phase formalism. We demonstrate the effectiveness of our method by performing 
comprehensive numerical tests on a model Ag/MgO/Ag heterostructure. 

PACS numbers: 71.15.-m 



I. INTRODUCTION 



The dielectric response of a metal-insulator interface 
to an applied electric field is a subject of great interest 
nowadays, as the ongoing miniaturization of electronic 
devices is reaching the atomic scale. In this regime, the 
properties of thin oxide films (used e.g. in nonvolatile 
ferroelectric memoriesi and as gate oxides in MOSFET 
transistors^) start to deviate from those predicted by 
macroscopic models, and cannot be disentangled from 
the metallic or semiconducting contact'^. Traditionally, 
deviations from bulk dielectric behavior in thin film ge- 
ometries have been described using phenomenological 
models, which have provided valuable insights into the 
qualitative consequences of low dimensionality and in- 
terfacing to a metallic electrode. However, such models 
do not have predictive quantitative capabilities for the 
many physical factors, such as interface chemistry^, fi- 
nite screening length in the metal^, or the presence of 
structural and point defects related to the growth pro- 
cess, that contribute to the observed properties. In order 
to separate these effects and to help improve over the 
existing models, the electronic and ionic response of a 
metal-insulator heterostructure to an external bias po- 
tential must be understood at the quantum mechanical 
level. 

First-principles calculations, especially within the 
framework of density functional theory, have proven to 
be a powerful tool in the fundamental understanding of 
metal-ceramic interfaces®. On the other hand, ab-initio 
studies of interfacial dielectric properties are still in their 
infancy, and only insulating heterostructures have been 
successfully simulated"' in a finite external field. Indeed, 
when one of the two materials is a metal (as required 
for the simulation of a realistic thin-film capacitor), the 
presence of partially filled bands at the Fermi level pre- 
vents the straightforward application of the Berry-phase 
technique^, on which most finite electric field calcula- 
tions have been based so far^'^. In order to gain further 
insight into the rich domain of phenomena occuring at 
metal-dielectric junctions there is the clear need for a 
methodology that allows one to work around this obsta- 
cle. 



In this work we demonstrate how this problem can 
be solved by using techniques and ideas borrowed from 
Wannier-function theory, which is an appealing alterna- 
tive to the traditional Berry-phase formalisniiS . In par- 
ticular, we show that the insulating nature of a capacitor 
under an applied bias potential is reflected at the mi- 
croscopic level in the confinement of the metal-induced 
gap states, whose propagation in the insulator is forbid- 
den; this property allows for a rigorous definition of the 
macroscopic polarization in real space. The coupling to 
an external field is then introduced through the usual def- 
inition of the electric enthalpy, in close analogy to purely 
insulating systems— i^. 

We demonstrate the effectiveness of our technique us- 
ing a model MgO/Ag(100) heterostructure; the simple 
properties of this well studied metal-ceramic interface 
make it an ideal test case for our method. In addition 
to providing a technical validation to our approach, our 
results indicate that MgO films retain an ideal classical 
behavior in the ultrathin limit, corroborating the picture 
that emerged from previous investigations^'^. 

This work is organized as follows. In Section II we de- 
fine the polarization, the electric enthalpy and the Hamil- 
tonian for our new finite-field approach, and show how it 
can be used to calculate the capacitance and local permit- 
tivity profiles of a general metal/insulator heterostruc- 
ture. In Section III we relate our scheme to other possi- 
ble approaches. In Section IV we present the numerical 
test on Ag/MgO in full detail, with particular focus on 
the new properties of the system that are made accessi- 
ble by our scheme (dynamical charges at the interface, 
electronic response, etc.). In Section V we briefiy dis- 
cuss the applicability of our method to range of physical 
problems, and present the summary and conclusions. 



II. METHOD 

To place our developments in context, we begin by 
reviewing the discrete Berry-phase formalism introduced 
by King-Smith and Vanderbilt^ (KSV) for a general insu- 
latingcrystal; we then adapt the method we developed in 
Ref. [13 to obtain a real-space expression for the polariza- 
tion. Next (Sec. Ill B|) we show how this procedure can be 
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readily extended to periodic metal/insulator heterostruc- 
tures. Finally, in Sec. Ill CI we introduce the couphng to 
an external field in the total energy and the Hamiltonian, 
which completes the derivation of our finite-field method. 
A discussion of the first-principles capacitance per unit 
area is provided in Sec. IIID( and links our method to 
the theory of local permittivity^'^, which was originally 
developed for purely insulating heterostructures. 



A. Insulating systems 

We start with a periodic insulator described by three 
real space lattice vectors R^; for convenience, we im- 
pose that Ri = (a, 0,0) is perpendicular to R2,3 (the 
latter two lie therefore on the yz plane). The Bril- 
louin zone for such a Bravais lattice is a prism in re- 
ciprocal space defined by the reciprocal space vectors 
Gi — (27r/a, 0, 0) and 62,3. We choose a discrete k- 
point sampling of the type k = jb + kj^, where the vec- 
tor b/parallel — (27r/L, 0,0) = Gi/iVy spans a regular 
one-dimensional mesh of dimension TVy = L/a and kj^ 
belongs to a set of N± special points in the perpendic- 
ular plane; the total number of independent fc-points is 
then A^xiV||. The x component of macroscopic polariza- 
tion (P) is commonly evaluated on this reciprocal space 
mesh by first defining a set of matrices: 



2kPnk+b|| 



(1) 



e-^^it'-knk). In 



where it is understood that |unk+G|| 
the KSV approach^ (P) is then calculated as a discrete 
Berry phase for each individual "stripe" of fc-points: 

(P)Ber.y(ki) = -^^^ImlndetM'^^+^b; (2) 

j 

the resulting values are finally averaged over the set of 
special points kj^ with corresponding weights : 



{P)b 



e-rry — / ^ 
kx 



(k±). 



(3) 



As an alternative, within single-particle theories, the 
macroscopic polarization can be equivalently related to 
the centers of the Wannicr functions^. In Ref. [H we 
have introduced a two-step procedure to obtain accurate 
values of the polarization by using a real-space technique 
based on Wannier functions: 

• First, the maximally localized Wannier functions^ 
are determined iteratively starting from a set of 
ground-state Bloch orbitals; 

• Then, the charge distribution of the Wannier func- 
tion is constructed in real space, and the center 
defined by a saw-tooth function. 

In the following, we modify this strategy to operate in one 
dimension instead of three, for several reasons: i) this 



way we can achieve an even more transparent analogy 
between the Wannier-function real-space method and the 
KSV Berry-phase formula; ii) an external field couples 
to one component of the polarization only; iii) it can be 
generalized in a very natural way to a system which is at 
the same time a two-dimensional conductor and a one- 
dimensional insulator (see Section fll Bp . 

When the first step above is restricted to one dimen- 
sion only, the resulting "hermaphrodite" Wannier func- 
tionsi^J^ are maximally localized along the x direction 
and Bloch-like in the orthogonal yz plane and are much 
easier to obtain compared to their three-dimensional 
counterparts. In fact, the unitary transformation that 
yields such a set of localized orbitals is constructed by 
means of the parallel-transport technique in a one-shot 
singular value decomposition, without the need for the 
iterative minimization of the spread functional^. The 
parallel transport technique consists^'* of enforcing the 
particular representation of the orbitals in which = 
^mn^n^ : where X^-^ are complex numbers of the form 
A}^^ = e"^" and is Hermitian. A definition of the 
Wannier centers can be written in terms of the phases 



27r 



(4) 



Interestingly, this definition of the centers of the parallel- 
transported hermaphrodite Wannier functions bears an 
exact relationship with the KSV discrete Berry-phase po- 
larization: 



{P)Berry{-k^) = ^Y.^n' 



(5) 



Once the approximate values, x}^-^, of the Wannier 
centers are available, the "refinement" in real space by 
means of saw-tooth functions (second step above) can 
be carried out in close analogy to Ref. First, the 
charge density associated with the hermaphrodite Wan- 
nier function lui^kj^) is explicitly constructed in a super- 
cell which is a multiple of the primitive one (of total 
length L along x) and integrated over the yz planes to 
yield a one-dimensional function of periodicity L: 



(6) 



p^(x)= dydzK-(r)|^ 



Then, the refined Wannier centers Xn are defined as: 



Xn^= I Xpl^{x)dx. 



(7) 



Finally, the expression for the polarization (P) follows: 



(8) 



The advantages of using the real-space refinement step 
described above are twofold: 
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FIG. 1: Schematic representation of the projected band struc- 
ture of a metal-insulator heterostructure in the absence of 
(top) and with (bottom) electric field. An effective bias po- 
tential AV between one slab and the repeated images is ob- 
tained. 

• The real-space centers x^-^ are identical to the 
Berry-phase centers in the limit of A^n — > oo, 
but for a discrete mesh they differ; in particular, 
x^^ converge exponentially to the asymptotic limit, 
while xJj-L have a much slower, polynomial conver- 
gence (as ^^(A^i"^))^^. Therefore in practice a rela- 
tively coarse fc-point mesh is sufficient for the cal- 
culation of P, with a significantly reduced compu- 
tational effort. 

• The real-space expression for the polariza- 
tion allows for a straightforward extension to 
metal/insulator heterostructures and the introduc- 
tion of the coupling to a finite external field, which 
is the main goal of this work and will be treated in 
depth in the next subsection. 



B. Metal/insulator heterostructures 

We remind the reader that, in a metal/insulator het- 
erostructure, the usual reciprocal-space expression for 
(P) based on the Berry phase is incompatible with the 
fractional orbital occupancies which stem from the metal- 
licity of the electrodes; however, we will show hereafter 
that the real-space expression of the polarization de- 
scribed above can be adapted to provide a rigorous de- 
scription also in this case. We start from the same as- 
sumptions and definitions as in the insulating case, ex- 
cept that we set A^n — 1, i.e. only the two-dimensional 
special-points average on kj^ is performed in the Brillouin 
zone sampling. This is by no means a limitation, since 
for a physically meaningful study of a polarized capac- 
itor a necessary condition is that there be no coupling 



between the electrodes across the dielectric film, i.e. the 
corresponding one-particle bands must be virtually dis- 
persionless along the x direction. 

Incidentally, we note that "hermaphrodite" WFs, 
which are the basic ingredient of the method, are ideally 
suited to representing the inherently delocalized nature 
of electrons in the two-dimensionally conducting plates, 
while providing a fast spatial decay (which is crucial to an 
accurate definition of (P)-~) along the field direction (x). 
Unfortunately, the "localization" transformation, upon 
which the real-space formulation is rooted, cannot be per- 
formed straightforwardly here because of the fractional 
occupancies at the Fermi level: in the finite-temperature 
functional the total energy is no longer invariant upon a 
unitary transformation of the occupied states^^. To work 
around this problem and achieve a localized representa- 
tion along a;, we divide the problem into three regimes, 
which (as above) are identified independently for every 
kj^ point of the surface Brillouin zone (a schematic rep- 
resentation is shown in Fig.[T]): 

• The conduction states with zero occupancy are dis- 
carded from the computation, since they do not 
contribute to the ground state energy or charge 
density. 

• The submanifold of fully occupied valence states 
(with energy eigenvalue e„ < Ep — 6m Fig. [T|) 
is localized separately by means of the parallel- 
transport algorithnii^iii; this operation leaves the 
total energy unaffected in complete analogy to the 
purely insulating case. 

• The partially occupied states in the region of the 
Fermi level (shown as "0 < / < 2" in Figure [T|) lie 
in the energy gap of the dielectric; as a result they 
are quantum mechanically confined in the metallic 
slab without taking any further action, since they 
cannot propagate through the insulator. 

The centers, x^^ , of the resulting set of localized or- 
bitals are then calculated, analogously to the insulating 
case, by using Eq.[7](for the metal-induced gap states we 
set Xn = and assume by convention that the metallic 
slab is approximately centered at x — 0). The definition 
of the polarization is identical to the insulating case, ex- 
cept that we must take the orbital occupation numbers 
f^^ explicitly into account: 

n 

The multi-valuedness in the definition of x^-^ (mod- 
ulo a lattice vector) is removed by imposing x^-^ g 
[—L/2, L/2], with the origin a; = corresponding to the 
center of the metallic slab and L the cell length^. Since 
the band energies change their value during structural 
relaxation, the whole procedure has to be repeated at 
the beginning of every iteration. This means that the 
number of states assigned to each energy window is not 
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fixed during the simulation, but is free to vary, following 
the evolution of orbital eigenvalues. 

We note that the exact boundaries of the energy win- 
dows, in particular the energy separating the orbitals 
that are considered metal-induced gap states and those 
that are localized by means of the parallel-transport pro- 
cedure, is to some extent arbitrary; this choice, however, 
has no influence on the result (provided that in the lower 
window all states are fully occupied and that the middle 
window does not overlap with the conduction or valence 
band of the insulator). The method is therefore well de- 
fined. 



C. Electric enthalpy and Hamiltonian 

The coupling of a periodic system to a finite external 
field S is described, in both the metallic and insulating 
cases, by the electric enthalpy^'^: 

^Eks-^£-{P), (10) 

where Eks is the Kohn-Sham total energy, Q. is the vol- 
ume of the supercell and {P) is the macroscopic polar- 
ization along the field axis; in a parallel-plate capacitor 
this is always perpendicular to the electrode plane. 

We write the contribution AH^ to the Hamiltonian 
due to the electric field operator as a non-local, state- 
dependent potential Fn{x) applied to each WF, Wn, in- 
dividually (we drop the superscript in the following 
discussion) : 

Ai?^ =f ^|F„w„)(w„|. (11) 

n 

The Fn{v) are periodic sawtooth functions^ of the 
form F„(r) = x, with x € [xn — L/2,Xn + L/2], i.e. 
-Fn(0,y, z) = for every n, and is discontinuous in a 
point located far from the center x„, where the corre- 
sponding orbital has vanishing importance. We note that 
the average value of F„{x) over the cell is irrelevant for 
a purely insulating system, but is crucial when a metal- 
lic electrode is present. It can be easily verified that 
our choice respects the correct limit in the case of an 
isolated metallic slab in vacuum, where the field opera- 
tor becomes local. This also means that, when a cen- 
ter Xn crosses the point x — L/2 (the center of the 
insulating slab), the corresponding polarizing potential 
£Fn{x) undergoes a discontinuous, rigid jump of mag- 
nitude AV = zL£L. For a sufficiently thick insulating 
layer, and as long as A is small enough not to bring va- 
lence or conduction states into the partially occupied re- 
gion close to the Fermi level, this jump leaves the system 
unaffected apart from a trivial shift of the eigenvalues 
corresponding to fully-occupied Bloch states (the discon- 
tinuity can also be effectively smoothened by smearing 
the Fn with Gaussian functions; a practical procedure to 
do this is explained in the Appendix). Interestingly, the 
large-field Zener instability reported for purely insulating 



bulk materials^'^ translates here to a Schottky-tunneling 
instability that has to be avoided by careful analysis of 
the band alignment. 

AH^ has in general a small anti-Hermitian part due to 
the fact that the Wannier functions are never exactly zero 
at the cell boundary but retain a small (exponentially 
vanishing) tail: 

{w,\AH^\wj) - {{w,\AH^\wi)y = {w,\{Fj-F,)\w,). 

The influence of this anti-Hermitian component is to in- 
duce a "unitary torque" within the manifold of the occu- 
pied orbitals, which is undesirable since, when working 
in the ensemble density functional theory framework (or 
even more clearly in the case of insulating systems) , the 
total energy (and the electric enthalpy) are both invari- 
ant with respect to such a unitary transformation; this 
is therefore a spurious contribution that we remove by 
introducing a correction term of the form: 

AH'\w,) = 1^ k,)K|(F, - F,)\w,). (12) 
j 

This term goes exponentially to zero with increasing in- 
sulator thickness, and vanishes if each Wannier-like state 
is strictly localized in a; € [xn — L/AjXn + L /4] . With this 
correction, we could always achieve optimal convergence 
to the electronic ground state. 

We want to stress here that, at the origin of the non- 
zero anti-Hermitian term fEg. I12p . there's an important 
physical feature, i.e. the tunneling current across a thin 
insulating film is never exactly zero (although it is expo- 
nentially decreasing as a function of film thickness). In 
order to study the metastable (long-lived) state of the 
polarized capacitor within a ground-state theory, we are 
obliged to make an approximation and discard the tun- 
neling current by truncating the tails of the evanescent 
conduction states induced by the metallic electrodes. We 
expect this approximation to be very good in cases where 
the tunneling current is so negligibly small that it does 
not contribute appreciably to the physics of the system. 
This condition is indeed satisfied in the regime where our 
method is exact, i.e. when the tails of the metal-induced 
gap states are confined to a region much smaller than 
the thickness of the insulating film, and therefore direct 
tunneling is suppressed. 

D. Local dielectric response and capacitance 

The theory of local permittivity that was originally de- 
veloped for purely insulating interfaces'"^ is directly ap- 
plicable to metal/insulator heterostructures, and can be 
readily generalized to introduce the first-principles ca- 
pacitance density per unit area, C. We start from the 
usual textbook definition of the capacitance density. 
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with AV^ = £L the potential drop between the electrodes, 
and a the surface density of the free charge stored in 
the device. From classical electrostatics, the free charge 
is equal to the value of the induced polarization deep 
inside the electrode, where the applied field is completely 
screened. In our geometry this occurs in the center {x — 
0) of the metallic slab, and Eq. [T3]can be rewritten as 



C 



AP{x = 0) 
£L ' 



(14) 



where AP{x) is the planar- and macroscopically- 
averaged^^ induced polarization. Finally, A_P(0) can be 
extracted from the relationship, derived in Ref. ll^, be- 
tween the local permittivity, €{x), and the induced polar- 
ization, and the corresponding globally averaged quanti- 
ties (e) and (AP>; 



1 X AF(x) 
{e)J (AP) 



(15) 



Since deep in the metal e{x) diverges [i.e. l/e(0) — 0], 
Eq. [15] reduces simply to AP{x = 0) = (AF) + £/(47r), 
and Eq. [13] reduces to 



C = 



(16) 



where (e) = An{AP)/£ 1 is the overall dielectric con- 
stant of the supercell, and (AP) is the polarization in- 
duced by the field £. We use equation [TBI to obtain the 
capacitance values in Section IIVI 

In order to appreciate the influence of the interface 
on the properties of the device, it is often interesting to 
compare the ab-initio value of the capacitance (Eq. [TB]) 
to the well-known classical formula, which is valid in the 
macroscopic regime and depends only on the thickness t 
and on the bulk permittivity of the dielectric film ei,uik'- 

Sbulk (17) 



C, 



Classical 



AlTt 



In a size regime comparable to the electronic wavelength 
the "thickness" of the insulating slab ceases to be a rigor- 
ous concept, since the transition from one bulk material 
to the other is not sharp, and their individual properties 
are strongly modified in close vicinity of the junction. 
In most cases we can define a "nominal" thickness as 
t — mto, where m is the number of layers of the crys- 
talline film and is the equilibrium interlayer spacing of 
the dielectric film far from the interface. Together with 
this definition of the thickness, we will use Eq.[T7]in Sec- 
tion [IV] to define a nominal capacitance Cq. 



III. ALTERNATIVE APPROACHES 

A. Metal/insulator heterostructures 

In principle, to address the problem of the dielectric 
transition at a metal/insulator interface, several alter- 
natives to our method can be explored; we will briefiy 
review in this section the most obvious solutions. 
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FIG. 2: Schematic representation of the self-consistent 
electrostatic potential (red curve) in unsupported 
metal/insulator/metal (a) and metal/insulator (b) slabs 
when an external field E^xt is switched on. In the geometry 
(a) by classical electrostatics the internal macroscopic field 
acting on the dielectric Sint is zero; in case (b) it corresponds 
to £ext/^buik, where etuik is the bulk dielectric constant of 
the insulator. 



We have already mentioned that the Berry phase for- 
malism is difficult to apply because of the partially oc- 
cupied states at the Fermi level. Conventional linear- 
response^ techniques seem also problematic, because 
the response to a uniform external field is again incom- 
patible with the overall metallic character of the sys- 
tem. A promising alternative is provided by quantum 
transport techniques, either within the noncquilibrium 
Green's function formalism or within time-dipendent 
density functional theory2£; our approach should yield 
comparable results in the limit of small tunneling cur- 
rent, although probably at a substantially lower compu- 
tational cost. 

A different approach consists in introducing a layer of 
vacuum in the heterostructurc, in order to allow for the 
straightforward application of a local saw-tooth poten- 
tial; since a capacitor is a finite object this seems also 
a natural choice for the geometry (Fig. [2^). As the two 
electrodes lie in the same supercell, however, there is only 
one Fermi level that defines the chemical potential of an 
electron in both metallic slabs; in other words, an ex- 
ternal perturbation can not induce a potential difference 
between the capacitor plates. Instead, the system be- 
haves as a short circuit, analogously to a thick, homoge- 
neous metallic slab, irrespective of the dielectric material 
sandwiched in between. In this case, the external field 
is completely screened by the accumulation/depletion of 
charge at the free surfaces of the conductor. 

A possibility to work around this problem is to simu- 
late just "half" of the capacitor, i.e. a metal/insulator 
slab isolated in vacuum (Fig. [2]d). Again, an external 
saw-tooth function can be used both to apply a uniform 
external field and to define the polarization (P) ; from the 
self-consistent electronic ground state in this perturbed 
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potential we can extract charge density differences and 
local polarization profiles, that can then be used to es- 
timate the capacitance and interfacial properties. There 
are, however, some limitations to this approach: 

• The termination of the insulating lattice must not 
be metallic or we recover case 1 above; this rules out 
the simulation of most polar interfaces [e.g. (100)- 
oriented 5-1 or 3-3 perovskites] . Furthermore, even 
for "insulating" crystal surfaces, sometimes surface 
states due to dangling bonds substantially reduce 
the electronic gap, thus limiting the range of elec- 
tric field values that can be used in practice. 

• By classical electrostatics, the electric field in the 
vacuum is larger than the electric field in the insu- 
lator by a factor of e. This means that, for high- 
permittivity dielectrics, simulations are restricted 
to extremely small values of the internal field, and 
therefore the non-linear regime is not practically 
accessible by means of this strategy. Even worse, 
when the insulator is a ferroelectric, it has already 
been shown that the presence of a free surface 
makes it difficult to impose the correct electrical 
boundary conditions inside the film^^. 

The advantage of our method is that its limitations 
are exclusively dictated by fundamental physics (Schot- 
tky or direct tunneling), not by artifacts that must be 
introduced ad hoc in an unconventional simulation cell. 
We note that the band alignment at the metal/dielectric 
junction plays a crucial role in determining the properties 
of the interface; particular care must be used to assess the 
reliability of the exchange and correlation functional in 
reproducing this feature, at least at a qualitative level. 



B. Pure insulators 

A real-space approach for applying a finite external 
field to a periodic system was first proposed by Nunes and 
Vanderbilt^^ and applied to a tight-binding Hamiltonian, 
which was based on the linear-scaling scheme of Mauri, 
Galli and Car-^'^. The first application of this approach to 
a self-consistent density-functional calculation was pre- 
sented by Fernandez, Dal Corso and A. Baldereschi^l, 
and validated by the study of the dielectric properties of 
bulk Si and GaAs. The great advantage provided by the 
linear scaling functional (Ref. [23l ) is that it leads natu- 
rally to a description of the electronic structure in terms 
of Wannier functions which are strictly confined in space; 
this means that the real-space position operator is for- 
mally well-defined. The price to pay, however, is that 
one has to deal with a non-conventional functional, where 
the orthonormality constraint is not exactly satisfied and 
that might pose convergence problems during the search 
for the electronic ground state; we are not aware of prac- 
tical follow-ups of these techniques. 

The vast majority of the finite-field calculations re- 
ported so far were performed within the scheme of Souza, 
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FIG. 3: High-frequency (red dashed curve) and static (black 
shaded curve) inverse permittivity profiles for the Ag/MgO 
heterostructure. Bulk values for the inverse permittivity are 
indicated as thin dotted lines. The approximate positions of 
the atomic planes are sketched in the upper part of the figure. 



ifiiguez and Vanderbilt- or Umari and Pasquarello^, both 
based on the Berry-phase definition of polarization. In 
the case of pure insulators, our scheme can be considered 
as an interesting practical alternative to the latter two 
approaches, since it can yield substantially better accu- 
racy for the same choice of computational parameters. 



IV. APPLICATION AND NUMERICAL TESTS 

We demonstrate the effectiveness of our scheme by pre- 
senting resuhs for a model [Ag(100)]„/[MgO(100)]m het- 
erostructure, with the thicknesses of the Ag and MgO 
slabs set to n = 9 and to = 13 atomic layers respectively. 
We choose to fix the in-plane lattice constant to our cal- 
culated equilibrium value of bulk MgO (ao=7.86 a.u.), 
and apply a 3% uniform in-plane expansion to the Ag lat- 
tice in order to obtain a pseudomorphic structure. The O 
atoms are placed on top of the interface Ag atoms, which 
is the most favorable alignment; the interface Mg atoms 
then sit above a four-fold hollow site of the Ag(lOO) sur- 
face. 



Computational parameters and zero-field 
equilibrium structure 



Our calculations were performed within the local den- 
sity approximation^^ and the PAW methodSS, as imple- 
mented in our "in-house" code, with a plane wave cutoff 
of 40 Ry, a Gaussian smearing of 0.2 eV for the orbital oc- 
cupancies and 10 special fc-points in the surface Brillouin 
zone. For all electronic and ionic relaxations we used an 
extension of the Car and Parrinellc^^ method to metallic 
systems^SiS^. The out-of-plane atomic positions and cell 
parameter were relaxed in zero field until the maximum 
force was lower than 0.05 meV/A, and the spacing to 
between the central MgO layers was within 0.2% of the 
bulk equilibrium value. 
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TABLE I: Callen dynamical charges at the MgO/Ag(100) in- 
terface; in MgO the first values refers to Mg, the second to 
O. 



FIG. 4: (i) X (longitudinal) component of the Callen dynam- 
ical charges in the Ag/MgO heterostructure in units of e. (ii) 
Atomic relaxations induced by a field of 10~^ a.u. 



B. High-frequency response 

We then imposed an external field of 51.4 V/fj,m and 
relaxed the electrons to the ground state while keeping 
the ions fixed. The resulting high-frequency permittivity 
profile is shown as a thick dashed curve in Fig. [31 On 
the metal side of the interface the inverse permittivity 
drops quickly to zero, as expected, while it is very close 
to our calculated bulk MgO value after only three atomic 
layers on the insulator side (the bulk calculations were 
performed on a tetragonal 8-layer non-primitive super- 
cell with the same in-plane lattice parameter and /c-point 
sampling as the interface system; the out-of-plane spac- 
ing was set to to)- The calculated value in the middle of 
the MgO slab is e°°(L/2)=3.095, which is within 0.03% of 
our calculated bulk value, e^if.=3.Q9A. This level of ac- 
curacy reflects the favorable convergence of the real-space 
expression of the polarization, which was demonstrated 
m Ref. m. As a comparison, in a previous calculation 
of insulating heterostructures based on the Berry-phase 
approach, errors introduced solely by the finite field for- 
malism were estimated to be of the order of 10%^^. The 
overall dielectric constant of the heterogeneous supercell 



was (e 



= 5.203. 



C. Dynamical charges 

The ionic forces // extracted from the finite-field 
ground state calculation provide the Born effective 
charges = fi/£ of all atoms at once^; their values, 
however, are not a local property of the interface as they 
depend on the overall geometry and composition of the 
supercell. It is physically more insightful, instead, to look 
at the Callen (or longitudinal) dynamical charges, which 
are obtained by dividing the Born charges by the average 
high-frequency dielectric constant of the heterostructure 



(e°°), and are indeed local^. The Callen charges for 
the present system are represented in the upper panel of 
Fig. m their values for the atoms closest to the interface 
are reported in Tab. HI It is apparent that, already three 
or four layers away from the interface, the charges asso- 
ciated to the Mg and O ions have already converged to 
the bulk value; their absolute value is slightly reduced 
in the first two MgO layers, while at the same time the 
interfacial Ag atom has a non-negligible positive charge. 
This behavior is consistent with the penetration of the 
metallic states in the insulator, which partially screens 
the ionic response of the ions closest to the interface; at 
the same time, the external field has a finite penetration 
length in the metallic slab, which affects the interface Ag 
layer. 



D. Static response 

The ionic contribution to the static polarization can 
be either evaluated in the linear (zero field) approxima- 
tion with standard linear-response technique&iS., or non- 
perturbatively, by performing a direct relaxation in a 
finite external field^; we tested both approaches here. 
First we computed the symmetry-restricted force matrix 
of the system by the finite-difference method in zero field 
using small displacements of 0.0004 A along x. The varia- 
tion of the local polarization upon atomic displacements, 
SP{x,{Ri})/SRi, was also simultaneously extracted for 
each degree of freedom and used, together with the force 
matrix and the Born effective charges, to construct the 
static permittivity profile, e^{x), which is shown as a 
thick solid curve in Fig. [3] The behavior of the static 
and high frequency permittivities are similar, the only 
relevant difference being their magnitudes in the middle 
of the insulating slab, with e"(L/2) converging to a value 
of 9.227 (our calculated e°^;j, = 9.229). In the lower panel 
of Fig.|4]we show the field-induced atomic displacements 
(the center of the metallic slab has been assumed as a 
reference). As expected, the metallic slab remains prac- 
tically unperturbed, while interestingly, with this choice 
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of origin, both Mg and O sublattices relax opposite to 
the field direction; the relaxation of the Mg ions is much 
less important than the relaxation of the oxygens. 

Next, we performed a direct relaxation of the atomic 
positions in an external bias Ay = 1.15 V, which cor- 
responds roughly to a self-consistent field of 400 V//im 
throughout the dielectric; the resulting average permit- 
tivity was within 0.05% of the zero field linear response 
value, indicating as expected the absence of anharmonic 
effects in this field regime within numerical accuracy. 

E. Capacitance density 

Having assessed the reliability of our technique, and 
obtained accurate profiles, we next extracted the static 
and high-frequency values of the capacitance, obtaining 
C° = 31.2 fF//im2 and C°° = 10.3 fF/^m^ respectively. 
These values are respectively 3.1% and 1.2% higher than 
the classical estimates based only on the knowledge of 
the bulk permittivity and of the thickness of the film: 

(0,oo) 

^(0,oo) _ ^bulk /ION 
classical ^^^^^ ( ) 

The close similarity between C and Cdassicai indicates 
that the influence of the metal on the dielectric proper- 
ties of the MgO film is very small. Previous investigations 
of ultrathin MgO films deposited on Ag have already evi- 
denced a surprisingly fast recovery of the bulk oxide local 
electronic structure as a function of film thickness^i; here 
we see that this result also holds for the dielectric prop- 
erties of MgO films, which show a nearly ideal classical 
behaviour down to the nanometer scale. 
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FIG. 5: (i) Energy conservation during molecular dynamics in 
an external bias of 1.15 V. Represented are (black solid), 
+ Kion (blue dashed) and E^ + Kion + Kfake (red solid). 
The zero of the energy scale is set at the ground state total 
energy in zero field, (ii) Polarization, (P), as a function of 
time, (iii) Frequency-resolved dipolar activity of the system 
(black shaded curve); zone-center transverse optical phonon 
frequency of bulk MgO is also shown (red dashed line). 



spectrum broadened by a Gaussian of width 20 cm"-'^ and 
weighted by the dipolar activity of the mode, [Z*(cl')]^, 
where Z*{uj) is the mode effective charge^. The domi- 
nant peak, centered at 400 cm~^, closely corresponds to 
the main oscillation period of the polarization (P) during 
the dynamics; the frequency of the calculated zone-center 
transverse optical mode of bulk MgO (420 cm"""^) is indi- 
cated as a thick dashed line. 



F. Molecular dynamics 

Since the formalism is not exactly variational, it is 
important to assess its suitability for calculating finite- 
temperature effects by monitoring energy conservation 
during a microcanonical molecular dynamics run^*. We 
start from the centrosymmetric zero-field equilibrium 
structure and let the atoms evolve along the x axis un- 
der an external bias of AV — 1.15 V. The evolution of 
the potential energy , of the "physical" constant of 
motion + Kion and of the approximate "mathemati- 
cal" constant of motion E^ + Kion + Kfakc {Kion is the 
ionic kinetic energy, while Kfake is the fictitious elec- 
tronic kinetic energy2I,) are shown in the upper panel of 
Fig. [5l The quality of the energy conservation is excel- 
lent; the curve has no appreciable drift, even when en- 
larging the vertical scale by three orders of magnitude. 
In the lower panel we show the evolution of (P), which 
undergoes quasiperiodic oscillations of approximate fre- 
quency 390 cm^^. To link the molecular dynamics run to 
the lattice-dynamical properties of the system extracted 
from the force matrix, in the inset we show the phonon 



V. CONCLUSIONS 

In summary, we have developed a first-principles 
method for calculating the response of hybrid 
metal/insulator structures to an electric field, and 
have demonstrated its applicability by calculating the 
local permittivity and capacitance of a model thin-film 
capacitor. The method is also appropriate for insulat- 
ing materials at finite bias, where it provides greater 
accuracy than conventional techniques, in both the 
linear and anharmonic regimes with respect to field 
and temperature; this is crucial in the simulation of 
ferroelectric devices, where anharmonic effects are an 
essential physical ingredient for a realistic description"^ 3.. 
The generality of the technique is evidenced by the 
high quality of overall energy conservation in molecular 
dynamics, which suggests fascinating applications to 
phenomena occuring at liquid/electrode interfaces, e.g. 
the recently observed freezing of interfacial water under 
electric fields"^^. 
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APPENDIX A: NUMERICAL 
IMPLEMENTATION 

In this Appendix we will describe the technical de- 
tails of the implementation of our method within the 
projector-augmented wave^S (PAW) formalism. The 
same prescriptions apply to the closely related ultra- 
soft pseudopotential technique^, while the case of norm- 
conserving pseudopotentials is easily recovered by remov- 
ing the augmentation terms in the following expressions. 

1. Notation 

The basic ingredients of the PAW method are a set 
of all-electron (AE) and pseudo (PS) partial waves, re- 
spectively \4>T,i) and \4>T,i), and the corresponding pro- 
jector functions \pT,i)- The index i here embeds three 
indices: two angular momentum variables I and m, and 
a third one, n, that accounts for the fact that several 
partial waves per angular momentum can be used; the 
dependence on the atomic species/site is included in r. 
It is useful to define two kind of integrals based on these 
quantities, namely the overlap matrix elements: 

Qt.ij = (0r,i|0rj) - {4'T,t\4>T,j), (Al) 

and the localized dipole moments of the augmented 
charges: 

dT,tj = {4>r,%\A4>r,i) - {4>r,%\A4>r,i)- (A2) 

Most of these integrals are zero because of the angular 
momentum selection rules. The projectors are used to 
construct the local density matrices: 

^■^.U ^ X! /«k(V'rik|Pi) (Pi l-iAnk) (A3) 
nk 

2. Polarization 

Next, following Ref. we first separate the total po- 
larization into three parts: 

P = P^on + + P""^. (A4) 

respectively the ionic contribution, the Berry phase (BP) 
contribution and the expectation value (EV) term. The 
first is trivial, and is given by the dipole moment of the 
bare pseudo-atomic charges: 



The third (EV) term is also easy to compute and is 
nothing but the total dipole moment of the augmented 
charges: 

P^"" =-^Y.'^r,,Wr^,,. (A6) 

The second (BP) term has the same form as Eq. [2] 

^P)llrry{^i:) = - 1 1- ^ Im In dct iW-- (A7) 

j 

however, the M matrix here must include the appropriate 
augmentation terms: 

^mn = (Wmk|Unk+b||) + 
E Qr,.,(^™k|p!''^)(pf Vnk+b,|) (A8) 

The parallel transport construction explained in 
Sec. Ill Al then provides the unitary rotations that localize 
the Wannier functions along the field axis, which lead eas- 
ily to the one-dimensional PS Wannier densities p„(x). 
Within PAW, however, there is one supplementary step 
that must be taken when explicitly constructing the local- 
ized orbitals in real space, i.e. one must reconstruct the 
localized augmented charges associated with the Wan- 
nier orbital; in particular we will need in the following 
the total augmented charge per site associated to the n- 
th Wannier function: 

9T,n = ^ QrAi, {Wn \Pi) iP] \Wn) (A9) 

3. Saw-tooth functions 

Finally, to complete our definition of the polarization 
we need to construct the saw-tooth functions Fn{x) in 
real space; in the following we present a recipe to do so 
while ensuring a continuous evolution of the Hamiltonian 
operator. Given a Wannier center Wn of approximate 
center Xn (we always set a;„ = for the metal-induced 
gap states with fractional orbital occupancies), defined 
in a one-dimensional cell of length L, we construct first 
the periodic function: 

oo 

Hr,{x) = 1 - A ^ g-[x-2„ + (i+l/2)L]V-= (AlO) 
/ — — oc 

where A = L/(y^(t) is such that the average value of H 
is zero and a is about 0.5-1 bohr. Then we define our 
saw-tooth function as the integral of i?„: 

Fn{x) = f Hn{y)dy (All) 
Jo 

We note that this definition ensures that Fn (0) = for all 
n; the finite value of a effectively avoids any discontinuity 
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in the evolution of the Hamiltonian while facilitating the 
plane-wave representation. Then the "refined" center x„ 
is defined as: 

Xn = / Fn{x)pn{x)dx + y2Fn{Xr)qr,n+Xn- Fn{Xn), 

Jo 

(A12) 

where q-r.n is the augmented charge at the atomic site 
T associated with the n-th Wannier function, and X-r is 



the atomic position. Note that the charges qr.n and the 
PS densities satisfy the usual generalized normalization 
condition: 

/ Pn{x)dx + y2qr.n ^ (A13) 

The Fn{x) functions are then used in Eq. [Tl]to build 
the electric field contribution to the Hamiltonian. 
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According to the modern theory of polarization, only dif- 
ferences of (P) between two states of the system con- 
nected by a continuous transformation are observable 
quantities, while the absolute value of (P) is only de- 
fined modulo a quantum of polarization^. The condition 
G [~L/2, L/2] (with the metallic slab approximately 
centered in the origin) , is important to ensure that all frac- 
tionally occupied state belong to the same supercell; in this 
case it is easy to show that (P) is well defined and has the 
same properties as in the insulating case. 
The Fn{r) are constructed in reciprocal space in order to 
avoid Fourier aliasing errors and preserve analytical gradi- 
ents. 



